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Abstract. A new description of the dynamics of warped accretion discs 
is presented. A theory of fully nonlinear, slowly varying bending waves 
^| is developed, involving a proper treatment of viscous fluid dynamics but 

. neglecting self-gravitation. The special case of resonant wave propagation 

found in inviscid Keplerian discs is not considered. The resulting descrip- 
tion goes beyond existing linear and nonlinear theories, and is suitable 
for numerical implementation. 
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1. Introduction 
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I Non-planar accretion discs are believed to exist in many astrophysical situa- 

tions. In some cases, such as the nuclear disc of the galaxy NGC 4258 or the 
circumstellar disc of (5 Pictoris, a warped profile is observed more or less directly. 
— ■ In others, including the X-ray binaries Her X-l and SS 433, tilting of the disc 

is proposed in order to explain certain observational characteristics. A more 
Q ■ complete discussion is given by Pringle (this volume). 

In accretion discs, unlike the case of galactic warps, the effects of self- 
gravitation may often be neglected. The first attempts to describe the time- 
dependence of a warped disc within the framework of viscous fluid dynamics 
(Petterson 1978; Hatchett, Begelman and Sarazin 1981) indicated that the warp 
would simply diffuse away on a viscous time scale. However, Papaloizou and 
%— i ' Pringle (1983) showed that these theories were oversimplified and failed to con- 

serve angular momentum as a result of serious internal inconsistencies. 

Special complications arise because of the existence of a resonance in thin 
discs that are both Keplerian (or very nearly so) and inviscid (or very nearly 
so). A distinction must therefore be made between the generic (non-resonant) 
case and the resonant case which occurs when the two conditions 
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<; H/r and a <; H/r (1) 



are both satisfied. Here Q is the angular velocity, k is the epicyclic frequency, 
H/r is the ratio of the semi-thickness of the disc to the radius, and a is the 
dimensionless viscosity parameter (Shakura and Sunyaev 1973). Whether the 
resonant case occurs in realistic astrophysical situations is a matter of some 
uncertainty. Current estimates suggest that the resonant case is not relevant to 
discs in X-ray binaries or active galactic nuclei, but may apply in protostellar 
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discs. However, the resonance is delicate and might be destroyed by the effects 
of self-gravitation, magnetic fields, or turbulence. 

The non-resonant case was examined by Papaloizou and Pringle (1983), who 
considered a Keplerian disc with a significant viscosity. They treated the warp 
as a slowly varying disturbance with azimuthal wave number m = 1, using linear 
Eulerian perturbation theory. The resulting equation for the warp is a complex 
linear diffusion equation, which has been used in applications (e.g. Kumar and 
Pringle 1985). 

Other cases were treated by Papaloizou and Lin (1994, 1995), who con- 
sidered an inviscid, but not necessarily Keplerian, disc. Again, linear Eulerian 
perturbation theory was used, and the authors assumed the warp to be a slowly 
varying normal mode of the disc. In the non-Keplerian (non-resonant) case the 
warp obeys a dispersive linear wave equation, while in the Keplerian (resonant) 
case the warp obeys a non-dispersive linear wave equation. By considering the 
effect of a small viscosity on the inviscid modes, Papaloizou and Lin connected 
this theory with that of Papaloizou and Pringle (1983), showing how the tran- 
sition occurs between wave-like and diffusive behaviour in Keplerian discs when 
a f» H/r. This theory has also been used in applications (e.g. Papaloizou and 
Terquem 1995). 

A difficulty arises with these theories because the Eulerian method is for- 
mally valid only when the tilt angle (3(r,t) of the warp satisfies \(5\ <C H/r, a 
condition violated - almost by definition - by any interesting, observable warp. 
Of course, a linear theory may be derived by any convenient method, and vi- 
olation of the above condition does not necessarily imply that the linearized 
dynamics ceases to apply, only that the Eulerian method cannot be used to 
justify it. Indeed, it is clear the Eulerian method is basically not well suited 
to the description of a warp or even a rigid tilt, and that some kind of La- 
grangian method ought to be used to describe warps of significant amplitude 
and to discuss nonlinear effects. The Eulerian method offers little information 
about nonlinear effects, although it does predict that, in the resonant case, an 
amplitude | d/3/d In r\ ~ H/r would result in horizontal shearing motions in the 
disc comparable to the sound speed, which are expected to be unstable (Kumar 
and Coleman 1993; Gammie, Goodman and Ogilvie, in preparation). 

Meanwhile, Pringle (1992) developed a different approach in which the 
forms of the equations governing a warped viscous disc are derived simply by 
requiring mass and angular momentum to be conserved, but without reference 
to the detailed internal fluid dynamics of the disc. In this scheme, neighbouring 
rings in the disc exchange angular momentum by means of viscous torques which 
are of two kinds. One kind of torque (associated with a kinematic viscosity co- 
efficient v\) acts on the differential rotation in the plane of the disc, and leads 
to accretion; the other kind (associated with v-i) acts to flatten the disc. The 
equations derived (originally by Papaloizou and Pringle 1983) are 

_ + -- (rEBr )=0 
for the surface density £(r,i), and 

+ i|p***) - It {^fe) + i| (i^a^l) (3) 
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for the angular momentum, in the absence of external torques. Here v r (r,t) 
is the mean radial velocity and £(r, t) is the tilt vector, which is a unit vector 
parallel to the orbital angular momentum of the ring. These equations constitute 
a straightforward generalization of the equations used for a flat disc (e.g. Pringle 
1981). 

This approach is ostensibly valid for warps of large amplitude and would 
therefore represent an advance on the linear theory. It has been used in several 
applications, in particular to identify the radiation-driven instability (Pringle 
1996) and to explore its linear theory (e.g. Maloney, Begelman and Pringle 
1996) and nonlinear evolution (e.g. Pringle 1997). However, because the equa- 
tions of Pringle (1992) are derived somewhat heuristically, without reference to 
the detailed internal fluid dynamics of the disc, some doubts remain over the 
validity of this method. For example, has any internal degree of freedom of the 
rings been neglected? Is the interaction between neighbouring rings purely of 
the assumed form of viscous torques? Are there any nonlinear fluid dynamical 
effects that might limit the amplitude of the warp? How are the viscosity coef- 
ficients v\ and U2 related? It is therefore important to investigate whether the 
equations of Pringle (1992) can be derived ab initio from the three-dimensional 
fluid dynamical equations, and to understand how they connect to the previously 
established linear theory. 

An entirely different approach to warps of large amplitude, based on three- 
dimensional numerical simulations, is described by Nelson (this volume). 

2. Outline of the method 

The basis of the approach is to develop a hydrodynamic theory of unforced 
bending waves of large amplitude in a thin disc. This is most naturally done in 
the case of a spherically symmetric external potential, in which a flat disc has 
no preferred plane of orientation. Continuous with the zero-frequency rigid-tilt 
mode, there exist bending waves - with azimuthal wave number m = 1 in linear 
theory - that are slowly varying in time and space. In contrast, most of the 
modes in a thin disc vary on a time-scale comparable to and on a length- 
scale comparable to H (e.g. Lubow and Pringle 1993). Therefore these bending 
waves are special and deserve a special analysis. 

A full account of the analysis has been given elsewhere (Ogilvie 1998). It will 
suffice here to explain the various steps involved with an emphasis on physical 
interpretation. 

2.1. Warped coordinates 

The limitations of an Eulerian approach to the fluid dynamics of a warped disc 
have already been noted. However, a fully Lagrangian approach is not suitable 
for dealing with a differentially rotating flow. Therefore a 'semi-Lagrangian' 
method is used, in which the equations are derived in Eulerian form but referred 
to a coordinate system that follows the principal warping motion of the disc. 
Warped spherical polar coordinates (r,6,4>) are defined as follows (Figure 1). 
r is the spherical radial coordinate. On each sphere r = constant the usual 
angular coordinates (6, <p) are defined, but with respect to an axis that is tilted 
to point in the direction of the unit vector £(r,t). The intention is that the 
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Figure 1. Warped spherical polar coordinates (r,9,4>). 

disc matter on each sphere will lie close to 9 = ir/2. This coordinate system is 
similar to that introduced by Petterson (1978) but has certain advantages. It 
is valid for any amplitude of warp, whereas Petterson's system, which is based 
on tilted cylinders, breaks down for large amplitudes or for insufficiently thin 
discs. The present system also has the same Jacobian determinant as ordinary 
spherical polar coordinates, which simplifies the equations to some extent. 

The equations of fluid dynamics can be readily derived in these coordinates 
using methods of elementary vector calculus. It is assumed that the fluid obeys 
the compressible Navier-Stokes equation with isotropic viscosity, but viscous 
heating is neglected. The principal features of the equations are (i) the ap- 
pearance of a number of fictitious forces and (ii) the introduction of a 'modified' 
radial derivative operator. The effect of the latter feature is that in a warped disc 
the vertical pressure gradient gives rise to a radial force which drives horizontal 
motions (Figure 2). 

2.2. Thin-disc asymptotics 

The assumption that the disc is thin allows progress to be made by introduc- 
ing scaled variables and asymptotic expansions. Although informal thin-disc 
approximations are common in accretion-disc theory, the rigorous procedures 
and wide applicability of asymptotic methods are not always appreciated. In 
the present problem, the basic small parameter e is a characteristic value of 
H/r. The problem is simplified by (i) separating the 'fast' orbital time-scale 
from the 'slow' time-scale of viscous and dynamical evolution of the warp, (ii) 
separating the 'fast' orbital velocity, the 'intermediate' oscillatory velocities in- 
duced by the warp and the 'slow' accretion velocity, and (iii) allowing thin-disc 
geometrical simplifications. For example, the relative velocity components have 
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Figure 2. The mechanism by which the stratification of pressure (and 
viscous stress) drives horizontal motions in a warped disc. 

the expansions 

v r (r,9,<l>,t) = ev rl {r,<t>,(,T) + e 2 v r2 (r,<j>,(,T) + 0(e 3 ), (4) 

v e {r,0,<l>,t) = ev ei (r,0,(,T) + e 2 v e2 (r,cj ) X,T) +0(e 3 ), (5) 

v^{r,e,<l>,t) = rO(r)sin^ + e^i(r,0,C,T)+e 2 ^ 2 (r,0,C,T)+O(e 3 ), (6) 

where ( = ^ 1 (ir/2 - 9) is a scaled dimensionless vertical coordinate in the 
disc, and T = e 2 t is the slow time coordinate. Note that all variables are non- 
axisymmetric except for the orbital angular velocity Q(r). 

2.3. Formal structure of the problem 

After substituting the asymptotic expansions into the dynamical equations, a 
number of equations are extracted at different orders. One of these determines 
the orbital angular velocity £l(r) as that of a free particle in a circular orbit 
in the given potential. The remaining equations may be divided into two sets. 
Physically, Set A determines the 'intermediate' velocities, while Set B determines 
the 'slow' velocities. Mathematically, Set A consists of coupled nonlinear partial 
differential equations (PDEs) in two dimensions (<f>,(), while Set B consists of 
coupled linear PDEs whose coefficients depend on the solution of Set A. 

It turns out that, while Set A must be solved in full, all the information 
required from Set B can be extracted by integration. Indeed, Set B can be 
manipulated to derive the basic conservation equations, 

» (jVsm) + 1 ^ y = id r QX7 3tf t \ + is ( Q2 i r w^) 

ot r or r or V / r or \ Or J 

+ l -^(Q 3 lrW£xf). (8) 
r Or \ Or J 
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for mass and angular momentum. Note that this angular momentum equation 
is similar to, but more general than, equation (|3|). Here I(r, t) is the azimuthally 
averaged second vertical moment of the density, which in cylindrical polar coor- 
dinates would be written 

and is an important dynamical quantity in the theory of bending waves. The 
quantities Qi are dimensionless coefficients which can be calculated from the 
solution of Set A, and which depend on the amplitude of the warp as well as 
the rotation law and details of the thermodynamics and viscosity. The angular 
momentum equation implies that each ring in the disc experiences torques of 
three kinds from its neighbours: coefficient Q\ represents a torque tending to 
spin up (or down) the ring. This would be the usual viscous torque proportional 
to dO/dr in a flat disc, but in a warped disc there is an additional contribution, 
not proportional to dO/dr, due to a correlation between the radial and azimuthal 
velocities induced by the warp; this vanishes in an inviscid disc because the radial 
and azimuthal velocities are perfectly out of phase. Coefficient Q2 represents a 
torque tending to align the ring with its neighbours, which acts to flatten the 
disc; this also vanishes in an inviscid disc. Coefficient Q3 represents a torque 
tending to make the ring precess if it is misaligned with its neighbours; this leads 
to the dispersive wave-like propagation of the warp. 

2.4. Separation of variables 

The further assumptions that (i) the disc is polytropic (or isothermal) locally in 
radius and (ii) the viscosity coefficients are locally proportional to the pressure 
allow Set A to be solved by separation of variables, since the vertical dependence 
of the solution can be determined by inspection. The problem is then reduced to 
solving a set of dimensionless nonlinear ordinary differential equations (ODEs) 
in azimuth, subject to periodic boundary conditions, and extracting the three 
dimensionless coefficients Qi, Q2 and Q3 from the solution. The system is 
of fourth order and involves five dimensionless parameters: the dimensionless 
amplitude of the warp, 

dl 
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dlnr 

the dimensionless squared epicyclic frequency, 

g- d y>, (id 

dlnr 

the polytropic coefficient T and the dimensionless viscosity coefficients a and 
(for shear and bulk viscosities, respectively). 

2.5. Evaluation of the coefficients 

Examples of the numerical evaluation of the coefficients are given by Ogilvie 
(1998). Of particular importance is the diffusion coefficient Q% this is shown 
in Figure 3 for the case of an isothermal, Keplerian disc without bulk viscosity 
(r = 1; k 2 = 1; a h = 0). It can be seen that the resonant behaviour associated 
with the limit a — > is diminished as the amplitude of the warp increases. 
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Figure 3. Contour plot of the coefficient Q2 for an isothermal, Kep- 
lerian disc without bulk viscosity. The horizontal and vertical coordi- 
nates are the dimensionless amplitude of the warp and the dimension- 
less viscosity parameter, respectively. 
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Figure 4. The three modes that couple to produce a cubic nonlinearity. 
2.6. Nature of the nonlinearity 

A weakly nonlinear theory can be developed analytically (Ogilvie 1998) which 
provides truncated Taylor series of the form a + b\ip\ 2 + 0(|V>| 4 ) for the three 
coefficients. The weakly nonlinear theory therefore contains a cubic nonlinearity 
which arises through a three-mode coupling. The modes involved (Figure 4) are 
(i) the 'tilt' mode or f mode (m = 1, of odd symmetry), consisting locally of a 
uniform vertical translation of the disc, (ii) an inertial or r mode {m = 1, of odd 
symmetry), consisting of a horizontal epicyclic motion proportional to z, and 
(hi) an acoustic or p mode (m = 2, of even symmetry), consisting of a vertical 
motion proportional to z. The compressive nature of the third mode explains 
why the nonlinear behaviour depends on the thermodynamics and, in principle, 
on the bulk viscosity. 

3. Discussion 

It has been shown that the nonlinear fluid dynamics of a warped accretion 
disc is described by conservation equations (Q) and (|8|) for mass and angular 
momentum, respectively. The two basic assumptions are that (i) the disc is thin 
and (ii) the resonant case [equation (|l])] does not occur. Equation (|8|) may be 
decomposed into an equation for the component of angular momentum parallel 
to £, 

,2( 



dr r Or 



d_i 2 

dr 



(12) 



and an equation for the tilt vector, 



Xr 2 n [ — +v T —)= Q{IrVL 2 — + -— Q 2 lr 3 fl 2 — 
^Ot Or ) or r or \ or J 



+ Q 2 lr 2 £l 2 
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e + -^[Q3lr 3 n 2 £x^) . (13) 
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The warp therefore satisfies a kind of wave equation of parabolic type, featur- 
ing advection, diffusion and dispersion, each with a nonlinear dependence on 
the amplitude of the warp. This form of equation arises for the following rea- 
son. According to Set A, the velocities induced by the warp are determined 
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instantaneously and locally in radius. These velocities are in a 'geostrophic' 
state in which inertial forces balance the pressure gradients and viscous forces; 
time-derivatives of the velocities in the inertial frame do not appear because 
of the assumed separation of time-scales. This means that the velocities do 
not constitute additional dynamical degrees of freedom. In the resonant case 
this approximation breaks down and the time-derivatives must be restored; the 
resulting equations for the warp are hyperbolic rather than parabolic. 

Since equation (||) is the angular momentum equation it is straightforward 
to add to the right-hand side any additional torque T due to tidal forcing, 
Lense-Thirring precession, radiation forces, self-gravitation, etc. 

Under the further assumptions of Section 2.4 concerning the thermodynam- 
ics and viscosity, the coefficients Qi can be consistently evaluated. For small- 
amplitude warps there is agreement with linear theory, but for larger amplitudes 
a numerical solution of the ODEs is required. 

The dispersion coefficient Q3 is generally smaller than the diffusion coeffi- 
cient Q2, and, to the extent that it can be neglected, the angular momentum 
equation (|3|) given by Pringle (1992) is valid. However, in general, both v\ and 
V2 depend on the amplitude of the warp are neither is equal to the usual verti- 
cally averaged viscosity v. The most important result from linear theory for a 
Keplerian disc (Papaloizou and Pringle 1983) was that, although v\ ~ u, 



Therefore, for typical estimates in the range 0.01 <^ a <J 0.1, a Keplerian disc 
is much more resistant to warping than would be estimated on the basis of 
the viscosity v. The nonlinear analysis shows that this resonant behaviour is 
diminished as the amplitude of the warp increases: the ratio z^/^i decreases with 
increasing amplitude of the warp, although it remains significantly larger than 
unity. The consequences of this and other complex properties of the system 
ought to be investigated using a numerical implementation of the governing 
equations incorporating a consistent determination of the three coefficients. 
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